R Notebook

Load frequently used packages

library(openxlsx)
library(ggrepel)
library(dplyr)
library(tidyr)
library(data.table)
library(broom)
library(broomExtra)
library(tibble)
library(sjstats)
library(car)
library(lme4)
library(lmerTest)
library(ggplot2)
library(tibble)
library(modelr)
library(tidyverse)
#library(miceadds)
library(ggforce)
require(openxlsx)
library(tidyverse)
library(caret)
library(glmnet)
library(ggplot2)
select <- dplyr::select
filter <- dplyr::filter

Print session info for reproducibility

save.sessionInfo <- sessionInfo()

save.sessionInfo

R version 4.1.1 (2021-08-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ComplexHeatmap_2.8.0 RColorBrewer_1.1-3   circlize_0.4.14      lmPerm_2.1.0         MASS_7.3-56          gridExtra_2.3       
 [7] glmnet_4.1-3         caret_6.0-91         lattice_0.20-45      ggforce_0.3.3        forcats_0.5.1        stringr_1.4.0       
[13] purrr_0.3.4          readr_2.1.2          tidyverse_1.3.1      modelr_0.1.8         lmerTest_3.1-3       lme4_1.1-29         
[19] Matrix_1.4-1         car_3.0-12           carData_3.0-5        sjstats_0.18.1       tibble_3.1.6         broomExtra_4.3.2    
[25] broom_0.7.12         tidyr_1.2.0          dplyr_1.0.8          ggrepel_0.9.1        ggplot2_3.3.5        openxlsx_4.2.5      
[31] skimr_2.1.3          data.table_1.14.2   

loaded via a namespace (and not attached):
  [1] readxl_1.4.0         backports_1.4.1      plyr_1.8.7           repr_1.1.4           splines_4.1.1        listenv_0.8.0       
  [7] TH.data_1.1-0        digest_0.6.29        foreach_1.5.2        htmltools_0.5.2      fansi_1.0.3          magrittr_2.0.3      
 [13] cluster_2.1.3        doParallel_1.0.17    tzdb_0.3.0           recipes_0.2.0        globals_0.14.0       gower_1.0.0         
 [19] matrixStats_0.61.0   sandwich_3.0-1       hardhat_0.2.0        colorspace_2.0-3     rvest_1.0.2          haven_2.4.3         
 [25] xfun_0.30            crayon_1.5.1         jsonlite_1.8.0       survival_3.3-1       zoo_1.8-9            iterators_1.0.14    
 [31] glue_1.6.2           polyclip_1.10-0      gtable_0.3.0         ipred_0.9-12         emmeans_1.7.3        GetoptLong_1.0.5    
 [37] sjmisc_2.8.9         future.apply_1.8.1   shape_1.4.6          BiocGenerics_0.38.0  abind_1.4-5          scales_1.2.0        
 [43] mvtnorm_1.1-3        DBI_1.1.2            Rcpp_1.0.8.3         xtable_1.8-4         performance_0.9.0    clue_0.3-60         
 [49] bit_4.0.4            stats4_4.1.1         lava_1.6.10          prodlim_2019.11.13   datawizard_0.4.0     httr_1.4.2          
 [55] ellipsis_0.3.2       pkgconfig_2.0.3      farver_2.1.0         nnet_7.3-17          sass_0.4.1           dbplyr_2.1.1        
 [61] utf8_1.2.2           tidyselect_1.1.2     rlang_1.0.2          reshape2_1.4.4       effectsize_0.6.0.1   munsell_0.5.0       
 [67] cellranger_1.1.0     tools_4.1.1          cli_3.2.0            generics_0.1.2       sjlabelled_1.1.8     evaluate_0.15       
 [73] fastmap_1.1.0        yaml_2.3.5           bit64_4.0.5          ModelMetrics_1.2.2.2 knitr_1.38           fs_1.5.2            
 [79] zip_2.2.0            future_1.24.0        nlme_3.1-157         xml2_1.3.3           compiler_4.1.1       rstudioapi_0.13     
 [85] png_0.1-7            reprex_2.0.1         tweenr_1.0.2         bslib_0.3.1          stringi_1.7.6        parameters_0.17.0   
 [91] nloptr_2.0.0         vctrs_0.4.0          pillar_1.7.0         lifecycle_1.0.1      GlobalOptions_0.1.2  jquerylib_0.1.4     
 [97] estimability_1.3     insight_0.17.0       R6_2.5.1             IRanges_2.26.0       parallelly_1.30.0    codetools_0.2-18    
[103] boot_1.3-28          assertthat_0.2.1     rjson_0.2.21         withr_2.5.0          S4Vectors_0.30.2     multcomp_1.4-18     
[109] bayestestR_0.11.5    parallel_4.1.1       hms_1.1.1            rpart_4.1.16         timeDate_3043.102    coda_0.19-4         
[115] class_7.3-20         minqa_1.2.4          rmarkdown_2.13       Cairo_1.5-15         pROC_1.18.0          numDeriv_2016.8-1.1 
[121] lubridate_1.8.0      base64enc_0.1-3     

Define in/out directory used in this script

dir <- "/Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data"
dir.Anno.GRSoriginal <- "/Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Annotation/GRSoriginal"
dir.Anno.GRSrevised <- "/Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Annotation/GRSoriginal"

Load MEGA to HTP ID key

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

MEGA.IDkey <- fread("MEGA_041822_MEGA2_to_HTP_ID_key_v0.1_JRS.csv")

MEGA.IDkey

Put the column names representing identifiers into a vector for reference below

ID_colnames <- c("FamilyID", "RecordID", "MEGA.FID", "MEGA.IID", "MEGA.LabID")

Load the analysis metadata

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

meta.MEGA.T21_visit1_Celiac <- fread("MEGA_041822_META_CeliacGRS_v0.1_JRS.csv")

meta.MEGA.T21_visit1_Celiac

For reference later, read in supplemental tables from Sharp et al., 2019

setwd(dir.Anno.GRSoriginal)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Annotation/GRSoriginal inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

Table_S2 <- read.xlsx("apt15826-sup-0001-supinfo.xlsx", sheet = "Table S2", startRow = 3) %>%
  separate(`Odds.Ratio.[95%.CI]`, into = c("OR", "CI"), sep = " ", extra = "merge", remove = TRUE) %>%
  rename(OR.GRS = OR,
         CI.GRS = CI,
         Weight.GRS = `Weight.(β)`)

Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1 rows [15].

Table_S2

setwd(dir.Anno.GRSoriginal)
Table_S3 <- read.xlsx("apt15826-sup-0001-supinfo.xlsx", sheet = "Table S3", startRow = 3) %>%
  rename(OR.GRS = OR,
         Weight.GRS = `Weight.(β)`)
Table_S3

Read in imputed HLA genotype data provided by Paul Norman

setwd(dir.GRSdata)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

hla.imputed <- openxlsx::read.xlsx("EspinosaHLAstudy_HLAcalls_PJN.xlsx")

#hla.imputed %>% colnames()
colnames(hla.imputed)<-c("MEGA.LabID",
                         "Race",
                         "Ethnicity",
                         "sample.id",
                         "HLA-A.1",
                         "HLA-A.2",
                         "HLA-B.1",
                         "HLA-B.2",
                         "HLA-C.1",
                         "HLA-C.2",
                         "HLA-DRB1.1",
                         "HLA-DRB1.2",
                         "HLA-DQA1.1",
                         "HLA-DQA1.2",
                         "HLA-DQB1.1",
                         "HLA-DQB1.2",
                         "HLA-DPB1.1",
                         "HLA-DPB1.2",
                         "Comments")

hla.imputed

# Check that MEGA.LabID == sample.id
hla.imputed %>% filter(MEGA.LabID!=sample.id | is.na(MEGA.LabID) | is.na(sample.id))

# Good

Join with the MEGA to HTP ID key and prepared analysis metadata

hla.imputed01 <- hla.imputed %>%
  # Remove sample.id since it is redundant, as shown above.
  select(-c(sample.id)) %>%
  # Remove metadata columns in the data we got back from Norman's group; to be replaced at time of analysis with latest version of metadata.
  select(-c(Race, Ethnicity)) %>%
  left_join(MEGA.IDkey, by = "MEGA.LabID") %>%
  select(colnames(MEGA.IDkey), everything()) %>%
  select(-c(HTP_participant, In_plink_fam)) %>%
  full_join(meta.MEGA.T21_visit1_Celiac, by = ID_colnames)

# View the data before further manipulation:
hla.imputed01

rm(hla.imputed); gc()

            used   (Mb) gc trigger   (Mb) limit (Mb)  max used   (Mb)
Ncells   5689714  303.9    8726178  466.1         NA   8726178  466.1
Vcells 162921203 1243.0  308334637 2352.5     204800 308334637 2352.5

print("Note: I am purposely calling the LabID column MEGA.LabID because it is an identifier linked to static data. The genetic data shouldn't vary by which sample was used for genotyping (if it does, that's most likely a quality control problem). In analyses using the MEGA genetic data, we will likely want to join the genetic datasets with metadata by the LabID relevant to the other key variables in the analysis. For example, in this project we are looking at Celiac status at visit 1, so when we add in our metadata we will probably want to use age at visit 1 instead of age at time of the sample used for genotyping.")

[1] "Note: I am purposely calling the LabID column MEGA.LabID because it is an identifier linked to static data. The genetic data shouldn't vary by which sample was used for genotyping (if it does, that's most likely a quality control problem). In analyses using the MEGA genetic data, we will likely want to join the genetic datasets with metadata by the LabID relevant to the other key variables in the analysis. For example, in this project we are looking at Celiac status at visit 1, so when we add in our metadata we will probably want to use age at visit 1 instead of age at time of the sample used for genotyping."

Convert HLA genotype dataset into long format

hla.imputed_long <- hla.imputed01 %>%
  gather(key="Locus", value="Allele", `HLA-A.1`:`HLA-DPB1.2`) %>%
  mutate(Locus = gsub("[.]1", "", Locus),
         Locus = gsub("[.]2", "", Locus)) %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID,
         Locus, Allele, Comments, everything())

hla.imputed_long

Verify that each participant has exactly two HLA alleles for each locus

if (hla.imputed_long %>%
  filter(EXCLUDE_from_analysis == 0) %>%
    select(MEGA.IID, Locus, Allele) %>%
    group_by(MEGA.IID, Locus) %>%
    summarise(N=n()) %>%
    arrange(desc(N)) %>%
    filter(N!=2) %>%
    nrow() == 0) {
  print("All IIDs have exactly 2 alleles for each locus. Good!")
  } else {
    print("Oops, not all IIDs have exactly 2 alleles for each locus (and they should).")
  }

`summarise()` has grouped output by 'MEGA.IID'. You can override using the `.groups` argument.
[1] "All IIDs have exactly 2 alleles for each locus. Good!"

if (hla.imputed_long %>%
  filter(EXCLUDE_from_analysis == 0) %>%
    select(RecordID, Locus, Allele) %>%
    group_by(RecordID, Locus) %>%
    summarise(N=n()) %>%
    arrange(desc(N)) %>%
    filter(N!=2) %>%
    nrow() == 0) {
  print("All RecordIDs have exactly 2 alleles for each locus. Good!")
  } else {
    print("Oops, not all RecordIDs have exactly 2 alleles for each locus (and they should).")
  }

`summarise()` has grouped output by 'RecordID'. You can override using the `.groups` argument.
[1] "All RecordIDs have exactly 2 alleles for each locus. Good!"

Double check that the EXCLUDE filter produces a dataset that includes only one set of HLA typing per RecordID

if ( hla.imputed_long %>%
       filter(EXCLUDE_from_analysis == 0) %>%
     select(RecordID, MEGA.LabID) %>%
     unique() %>%
     group_by(RecordID) %>%
     summarise(N_LabIDs = n()) %>%
     filter(N_LabIDs > 1) %>%
     nrow() == 0 ) {
  print("All RecordIDs have only one LabID. Good!")
  } else {
    print("Oops, not all RecordIDs have only one LabID.")
    }

[1] "All RecordIDs have only one LabID. Good!"

View the long format dataset so far

hla.imputed_long

Make a dataset of the HLA-DQ alleles and output to CSV and TSV

analysisData.alleles <- hla.imputed_long %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus, Allele, EXCLUDE_from_analysis) %>%
  mutate(Locus_Allele = paste(Locus, Allele, sep = "")) %>%
  select(-c(Allele)) %>%
  mutate(Has.Locus_Allele = 1) %>%
  group_by(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, EXCLUDE_from_analysis, Locus_Allele) %>%
  summarise(Dosage.Locus_Allele = sum(Has.Locus_Allele)) %>%
  ungroup() %>%
  arrange(desc(Dosage.Locus_Allele)) %>%
  spread(key = Locus_Allele, value = Dosage.Locus_Allele) %>%
  gather(key = "Locus_Allele", value = "Dosage.Locus_Allele", `HLA-A*01:01`:`HLA-DRB1*16:02`) %>%
  mutate(Dosage.Locus_Allele = ifelse(is.na(Dosage.Locus_Allele), 0, Dosage.Locus_Allele)) %>%
  separate(Locus_Allele, into = c("Locus", "Allele"), sep = "[*]", extra = "merge", remove = FALSE) %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus, Locus_Allele, Dosage.Locus_Allele, EXCLUDE_from_analysis, everything())

`summarise()` has grouped output by 'FamilyID', 'RecordID', 'MEGA.FID', 'MEGA.IID', 'MEGA.LabID', 'EXCLUDE_from_analysis'. You can override using
the `.groups` argument.

analysisData.alleles

if (analysisData.alleles %>%
  filter(EXCLUDE_from_analysis == 0) %>%
    select(MEGA.IID, Locus, Locus_Allele, Dosage.Locus_Allele) %>%
    group_by(MEGA.IID, Locus) %>%
    summarise(Total_Dosage.Locus = sum(Dosage.Locus_Allele)) %>%
    arrange(desc(Total_Dosage.Locus)) %>%
    filter(Total_Dosage.Locus != 2) %>%
    nrow() == 0) {
  print("All IIDs have exactly 2 alleles for each locus. Good!")
  } else {
    print("Oops, not all IIDs have exactly 2 alleles for each locus (and they should).")
  }

`summarise()` has grouped output by 'MEGA.IID'. You can override using the `.groups` argument.
[1] "All IIDs have exactly 2 alleles for each locus. Good!"

if (analysisData.alleles %>%
  filter(EXCLUDE_from_analysis == 0) %>%
    select(RecordID, Locus, Locus_Allele, Dosage.Locus_Allele) %>%
    group_by(RecordID, Locus) %>%
    summarise(Total_Dosage.Locus = sum(Dosage.Locus_Allele)) %>%
    arrange(desc(Total_Dosage.Locus)) %>%
    filter(Total_Dosage.Locus != 2) %>%
    nrow() == 0) {
  print("All RecordIDs have exactly 2 alleles for each locus. Good!")
  } else {
    print("Oops, not all RecordIDs have exactly 2 alleles for each locus (and they should).")
  }

`summarise()` has grouped output by 'RecordID'. You can override using the `.groups` argument.
[1] "All RecordIDs have exactly 2 alleles for each locus. Good!"

analysisData.alleles %>%
  split(., .$Locus_Allele) %>%
  lapply(dim) %>%
  unique()

[[1]]
[1] 221  10

  
analysisData.alleles

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

fwrite(analysisData.alleles, "MEGA_041822_Imputed_HLA_Alleles_v0.1_JRS.csv")
fwrite(analysisData.alleles, "MEGA_041822_Imputed_HLA_Alleles_v0.1_JRS.tsv", sep = "\t")

Make a dataset of the HLA-DQ allele groups and output to CSV and TSV

analysisData.alleleGroups <- hla.imputed_long %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus, Allele, EXCLUDE_from_analysis) %>%
  separate(Allele, into = c("AlleleGroup", "Protein"), sep = "[:]", extra = "merge", remove = FALSE) %>%
  mutate(Locus_AlleleGroup = paste(Locus, AlleleGroup, sep = "")) %>%
  mutate(Has.Locus_AlleleGroup = 1) %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus_AlleleGroup, Has.Locus_AlleleGroup, EXCLUDE_from_analysis) %>%
  group_by(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, EXCLUDE_from_analysis, Locus_AlleleGroup ) %>%
  summarise(Dosage.Locus_AlleleGroup = sum(Has.Locus_AlleleGroup)) %>%
  ungroup() %>%
  arrange(desc(Dosage.Locus_AlleleGroup)) %>%
  arrange(Dosage.Locus_AlleleGroup) %>%
  spread(key = Locus_AlleleGroup, value = Dosage.Locus_AlleleGroup) %>%
  gather(key = "Locus_AlleleGroup", value = "Dosage.Locus_AlleleGroup", `HLA-A*01`:`HLA-DRB1*16`) %>%
  mutate(Dosage.Locus_AlleleGroup = ifelse(is.na(Dosage.Locus_AlleleGroup), 0, Dosage.Locus_AlleleGroup)) %>%
  separate(Locus_AlleleGroup, into = c("Locus", "AlleleGroup"), sep = "[*]", extra = "merge", remove = FALSE) %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus, Locus_AlleleGroup, Dosage.Locus_AlleleGroup, EXCLUDE_from_analysis, everything())

`summarise()` has grouped output by 'FamilyID', 'RecordID', 'MEGA.FID', 'MEGA.IID', 'MEGA.LabID', 'EXCLUDE_from_analysis'. You can override using
the `.groups` argument.

# analysisData.alleleGroups <- hla.imputed_long %>%
#   select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus, Allele, EXCLUDE_from_analysis) %>%
#   separate(Allele, into = c("AlleleGroup", "Protein"), sep = "[:]", extra = "merge", remove = FALSE) %>%
#   mutate(Locus_AlleleGroup = paste(Locus, AlleleGroup, sep = "")) %>%
#   select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus_AlleleGroup, EXCLUDE_from_analysis) %>%
#   mutate(Has.Locus_AlleleGroup = 1) %>%
#   group_by(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, EXCLUDE_from_analysis, Locus_AlleleGroup) %>%
#   summarise(Dosage.Locus_AlleleGroup = sum(Has.Locus_AlleleGroup)) %>%
#   arrange(desc(Dosage.Locus_AlleleGroup)) %>%
#   spread(key = Locus_AlleleGroup, value = Dosage.Locus_AlleleGroup) %>%
#   gather(key = "Locus_AlleleGroup", value = "Dosage.Locus_AlleleGroup", `HLA-A*01`:`HLA-DRB1*16`) %>%
#   mutate(Dosage.Locus_AlleleGroup = ifelse(is.na(Dosage.Locus_AlleleGroup), 0, Dosage.Locus_AlleleGroup))

# analysisData.alleleGroups <- hla.imputed_long %>%
#   select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Locus, Allele) %>%
#   unique() %>%
#   separate(Allele, into = c("AlleleGroup", "Protein"), sep = "[:]", extra = "merge", remove = FALSE) %>%
#   mutate(Locus_AlleleGroup = paste(Locus, AlleleGroup, sep = ""),
#          dummy = 1) %>%
#   spread(key = Locus_AlleleGroup, value = dummy) %>%
#   gather(key = "Locus_AlleleGroup", value = "Dosage", `HLA-A*01`:`HLA-DRB1*16`) %>%
#   mutate(Dosage = ifelse(is.na(Dosage), 0, Dosage))

if (analysisData.alleleGroups %>%
  filter(EXCLUDE_from_analysis == 0) %>%
    select(MEGA.IID, Locus, Locus_AlleleGroup, Dosage.Locus_AlleleGroup) %>%
    group_by(MEGA.IID, Locus) %>%
    summarise(Total_Dosage.Locus = sum(Dosage.Locus_AlleleGroup)) %>%
    arrange(desc(Total_Dosage.Locus)) %>%
    filter(Total_Dosage.Locus != 2) %>%
    nrow() == 0) {
  print("All IIDs have exactly 2 alleles for each locus. Good!")
  } else {
    print("Oops, not all IIDs have exactly 2 alleles for each locus (and they should).")
  }

`summarise()` has grouped output by 'MEGA.IID'. You can override using the `.groups` argument.
[1] "All IIDs have exactly 2 alleles for each locus. Good!"

if (analysisData.alleleGroups %>%
  filter(EXCLUDE_from_analysis == 0) %>%
    select(RecordID, Locus, Locus_AlleleGroup, Dosage.Locus_AlleleGroup) %>%
    group_by(RecordID, Locus) %>%
    summarise(Total_Dosage.Locus = sum(Dosage.Locus_AlleleGroup)) %>%
    arrange(desc(Total_Dosage.Locus)) %>%
    filter(Total_Dosage.Locus != 2) %>%
    nrow() == 0) {
  print("All RecordIDs have exactly 2 alleles for each locus. Good!")
  } else {
    print("Oops, not all RecordIDs have exactly 2 alleles for each locus (and they should).")
  }

`summarise()` has grouped output by 'RecordID'. You can override using the `.groups` argument.
[1] "All RecordIDs have exactly 2 alleles for each locus. Good!"

analysisData.alleleGroups %>%
  split(., .$Locus_AlleleGroup) %>%
  lapply(dim) %>%
  unique()

[[1]]
[1] 221  10

analysisData.alleleGroups

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

fwrite(analysisData.alleleGroups, "MEGA_041822_Imputed_HLA_AlleleGroups_v0.1_JRS.csv")
fwrite(analysisData.alleleGroups, "MEGA_041822_Imputed_HLA_AlleleGroups_v0.1_JRS.tsv", sep = "\t")

Derive dosage of Celiac-associated heterodimers

Pietzak, Michelle M., et al. “Stratifying Risk for Celiac Disease in a Large At-Risk United States Population by Using HLA Alleles.” Clinical Gastroenterology and Hepatology: The Official Clinical Practice Journal of the American Gastroenterological Association, vol. 7, no. 9, Sept. 2009, pp. 966–71. PubMed, https://doi.org/10.1016/j.cgh.2009.05.028. “HLA-DQ2.5 = DQA105-DQB10201 (DR3), HLA-DQ2.2 = DQA10201-DQB10202 (DR7), HLA-DR5 = DQA105-DQB103 , HLA DQ8 = DQA103-DQB10302 .” (https://pubmed.ncbi.nlm.nih.gov/19500688/)

temp <- hla.imputed_long %>%
  select(ID_colnames,
         EXCLUDE_from_analysis, EXCLUDE_reason,
         Locus, Allele) %>%
  filter(Locus=="HLA-DQA1" | Locus=="HLA-DQB1") %>%
  group_by(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID,
           EXCLUDE_from_analysis, EXCLUDE_reason,
           Locus, Allele) %>%
  summarise(Dosage = n()) %>%
  ungroup() %>%
  arrange(Dosage) %>%
  mutate(DQA1_05 = ifelse(Locus=="HLA-DQA1" & grepl("05[:]", Allele)==TRUE, Dosage, 0),
         DQA1_03 = ifelse(Locus=="HLA-DQA1" & grepl("03[:]", Allele)==TRUE, Dosage, 0),
         DQA1_02 = ifelse(Locus=="HLA-DQA1" & grepl("02[:]", Allele)==TRUE, Dosage, 0),
         DQA1_0201 = ifelse(Locus=="HLA-DQA1" & grepl("02[:]01", Allele)==TRUE, Dosage, 0),
         DQB1_02 = ifelse(Locus=="HLA-DQB1" & grepl("02[:]", Allele)==TRUE, Dosage, 0),
         DQB1_0202 = ifelse(Locus=="HLA-DQB1" & grepl("02[:]02", Allele)==TRUE, Dosage, 0),
         DQB1_0302 = ifelse(Locus=="HLA-DQB1" & grepl("03[:]02", Allele)==TRUE, Dosage, 0),
         # DQ7.5= DQA1*05-DQB1*03. Does not matter what the second number (i.e 0501, 0505 both count) (https://pubmed.ncbi.nlm.nih.gov/19500688/)
         DQB1_03 = ifelse(Locus == "HLA-DQB1" & grepl("03[:]", Allele)==TRUE, Dosage, 0),
         DQB1_0301 = ifelse(Locus == "HLA-DQB1" & grepl("03[:]01", Allele)==TRUE, Dosage, 0))

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(ID_colnames)` instead of `ID_colnames` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.
`summarise()` has grouped output by 'FamilyID', 'RecordID', 'MEGA.FID', 'MEGA.IID', 'MEGA.LabID', 'EXCLUDE_from_analysis', 'EXCLUDE_reason',
'Locus'. You can override using the `.groups` argument.

temp %>%
  select(MEGA.LabID,
         EXCLUDE_from_analysis, EXCLUDE_reason,
         Locus, Allele, Dosage,
         DQA1_05, DQA1_03, DQA1_0201,
         DQB1_02, DQB1_0202, DQB1_0301, DQB1_0302) %>%
  unique()

temp2 <- temp %>%
  select(-c(Locus, Allele, Dosage)) %>%
  group_by(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, EXCLUDE_from_analysis, EXCLUDE_reason) %>%
  summarise(DQA1_05 = sum(DQA1_05),
            DQA1_03 = sum(DQA1_03),
            DQA1_02 = sum(DQA1_02),
            DQA1_0201 = sum(DQA1_0201),
            DQB1_02 = sum(DQB1_02),
            DQB1_0202 = sum(DQB1_0202),
            DQB1_0302 = sum(DQB1_0302),
            DQB1_03 = sum(DQB1_03),
            DQB1_0301 = sum(DQB1_0301)) %>%
  ungroup() %>%
  mutate( DQ2.5 = ifelse(DQA1_05>=1 & DQB1_02>=1, 1, 0),
          DQ8 = ifelse(DQA1_03>=1 & DQB1_0302>=1, 1, 0),
          #DQ2.2.archive040521 = ifelse(DQA1_02>=1 & DQB1_02>=1, 1, 0),
          DQ2.2 = ifelse(DQA1_0201 >=1 & DQB1_0202 >=1, 1, 0),
          DQ7.5 = ifelse(DQA1_05>=1 & DQB1_0301>=1, 1, 0) # SharpEtAl DQ7.5 (DQA1*05, DQB1*03:01)
  )

`summarise()` has grouped output by 'FamilyID', 'RecordID', 'MEGA.FID', 'MEGA.IID', 'MEGA.LabID', 'EXCLUDE_from_analysis'. You can override using
the `.groups` argument.

          # 7.5= DQA1*05-DQB1*03. Marisa Stahl: "Does not matter what the second number (i.e 0501, 0505 both count)" (https://pubmed.ncbi.nlm.nih.gov/19500688/)
          # Pietzak: "HLA-DQ2.5  DQA1*05-DQB1*0201 (DR3), HLA-DQ2.2  DQA1*0201-DQB1*0202 (DR7), HLA-DR5  DQA1*05-DQB1*03, HLA DQ8 DQA1*03-DQB1*0302."
          # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4297300/
          # "The very few patients, who are not DQ2.5, DQ2.2, or DQ8, are almost all DQ7.5 (DQA1*05, DQB1*03:01) (Karell et al. 2003)."
          # Karell K, Louka AS, Moodie SJ, Ascher H, Clot F, Greco L, Ciclitira PJ, Sollid LM, Partanen J. HLA types in celiac disease patients not carrying the DQA1*05-DQB1*02 (DQ2) heterodimer: results from the European genetics cluster on celiac disease. Hum Immunol. 2003;64:469–477. doi: 10.1016/S0198-8859(03)00027-2.
          #DQ7.5_PietzakEtAl = ifelse(DQA1_05>=1 & DQB1_0301>=1, 1, 0),

df.heterodimer_dosage <- temp2 %>% arrange(MEGA.LabID)

print("Note: HTP0073A2 and HTP0374A carried DQ2.2 in 5/2021 version of data cleaning.")

[1] "Note: HTP0073A2 and HTP0374A carried DQ2.2 in 5/2021 version of data cleaning."

# View the dataframe:
df.heterodimer_dosage

Summarize the data in its interim state:

df.heterodimer_dosage %>%
  group_by(DQ2.5, DQ8, DQ2.2, DQ7.5) %>%
  summarise(N = n())

`summarise()` has grouped output by 'DQ2.5', 'DQ8', 'DQ2.2'. You can override using the `.groups` argument.

print("Per Marsia Stahl, \"DQ7.5 with DQ2.2 is equivalent to DQ2.5 'in trans'\"")

[1] "Per Marsia Stahl, \"DQ7.5 with DQ2.2 is equivalent to DQ2.5 'in trans'\""

print("Observation: All participants who have DQ2.2+DQ7.5 (DQ2.5 in trans) also have imputed genotypes for DQ2.5.")

[1] "Observation: All participants who have DQ2.2+DQ7.5 (DQ2.5 in trans) also have imputed genotypes for DQ2.5."

print("Decision: For these participants with DQ2.2+DQ7.5+DQ2.2, consider them DQ2.5 in trans.")

[1] "Decision: For these participants with DQ2.2+DQ7.5+DQ2.2, consider them DQ2.5 in trans."

print("Therefore:")

[1] "Therefore:"

print("If DQ2.5 == 1 & DQ2.2 == 0 & DQ7.5 == 0, 'DQ2.5'")

[1] "If DQ2.5 == 1 & DQ2.2 == 0 & DQ7.5 == 0, 'DQ2.5'"

print("If DQ2.5 == 1 & DQ2.2 == 1 & DQ7.5 == 1, 'DQ2.5 or DQ2.5 in trans'")

[1] "If DQ2.5 == 1 & DQ2.2 == 1 & DQ7.5 == 1, 'DQ2.5 or DQ2.5 in trans'"

df.heterodimer_dosage %>%
  mutate(DQ2.5_notInTrans = ifelse(DQ2.5 == 1 & DQ2.2 == 0 & DQ7.5 == 0, 1, 0),
         DQ2.5_or_DQ2.5inTrans = ifelse(DQ2.5 == 1 & DQ2.2 == 1 & DQ7.5 == 1, 1, 0)) %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID,
         EXCLUDE_from_analysis, EXCLUDE_reason,
         DQ2.5_notInTrans, DQ2.5_or_DQ2.5inTrans) %>%
  unique() %>%
  group_by(DQ2.5_notInTrans, DQ2.5_or_DQ2.5inTrans) %>%
  summarise(N = n())

`summarise()` has grouped output by 'DQ2.5_notInTrans'. You can override using the `.groups` argument.

print("We have 9 samples that are marked as only DQ2.5 but not possibly DQ2.5 in trans. Meanwhile we have 27 that are DQ2.5 in trans but not DQ2.5. Good, this is how the code was supposed to work.")

[1] "We have 9 samples that are marked as only DQ2.5 but not possibly DQ2.5 in trans. Meanwhile we have 27 that are DQ2.5 in trans but not DQ2.5. Good, this is how the code was supposed to work."

Create a variable that reflects the combination of heterodimers carried by each participant

df.Heterodimer_Combo <- df.heterodimer_dosage %>%
  mutate(DQ2.5 = gsub("1", "DQ2.5", DQ2.5),
         DQ2.5 = gsub("0", "", DQ2.5)) %>%
  mutate(DQ8 = gsub("1", "DQ8", DQ8),
         DQ8 = gsub("0", "", DQ8)) %>%
  mutate(DQ7.5 = gsub("1", "DQ7.5", DQ7.5),
         DQ7.5 = gsub("0", "", DQ7.5)) %>%
  mutate(DQ2.2 = gsub("1", "DQ2.2", DQ2.2),
         DQ2.2 = gsub("0", "", DQ2.2)) %>%
  mutate(Heterodimer_Combo = paste(DQ2.5, DQ8, DQ2.2, DQ7.5, sep = "+")) %>%
  mutate(Heterodimer_Combo = gsub("[+][+]", "+", Heterodimer_Combo)) %>%
  mutate(Heterodimer_Combo = gsub("[+][+]", "+", Heterodimer_Combo)) %>%
  select(ID_colnames, EXCLUDE_from_analysis, EXCLUDE_reason,
         Heterodimer_Combo) %>% #, DQ2.5, DQ8, DQ2.2, DQ7.5) %>%
  unique() %>%
  mutate(Heterodimer_Combo = ifelse(Heterodimer_Combo == "+DQ8+", "DQ8", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "+DQ7.5", "DQ7.5", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "DQ2.5+", "DQ2.5", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "+", "", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "+DQ8+DQ2.2+", "DQ8 + DQ2.2", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "+DQ2.2+", "DQ2.2", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "DQ2.5+DQ2.2+", "DQ2.5 + DQ2.2", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "+DQ8+DQ7.5", "DQ8 + DQ7.5", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "DQ2.5+DQ8+", "DQ2.5 + DQ8", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "DQ2.5+DQ2.2+DQ7.5", "DQ2.5 + DQ2.2 + DQ7.5", Heterodimer_Combo),
         Heterodimer_Combo = ifelse(Heterodimer_Combo == "DQ2.5+DQ7.5", "DQ2.5 + DQ7.5", Heterodimer_Combo)) %>%
  left_join(df.heterodimer_dosage, by = c(ID_colnames, "EXCLUDE_from_analysis", "EXCLUDE_reason")) %>%
  #filter(Heterodimer_Combo == "") %>%
  mutate(Heterodimer_Combo = ifelse(Heterodimer_Combo == "", "X + X", Heterodimer_Combo)) %>%
  mutate(Heterodimer_Combo = ifelse(Heterodimer_Combo == "DQ7.5", "DQ7.5 + X",
                                    ifelse(Heterodimer_Combo == "DQ2.5", "DQ2.5 + X",
                                           ifelse(Heterodimer_Combo == "DQ2.2", "DQ2.2 + X",
                                                  ifelse(Heterodimer_Combo == "DQ8", "DQ8 + X",
                                                         Heterodimer_Combo)))))
df.Heterodimer_Combo

df.Heterodimer_Combo %>%
  select(Heterodimer_Combo) %>%
  unique() %>%
  arrange(Heterodimer_Combo)

# DQ2.2 + X             
# DQ2.5 + DQ2.2             
# DQ2.5 + DQ2.2 + DQ7.5             
# DQ2.5 + DQ7.5             
# DQ2.5 + DQ8               
# DQ2.5 + X             
# DQ7.5 + X             
# DQ8 + DQ2.2               
# DQ8 + DQ7.5               
# DQ8 + X
# X + X

Create a variable that represents HLA-DQ Genotype the same way as in the Sharp et al. paper

# Notes:
# Of the 15 HLA-DQ genotypes represented in the Sharp paper, only 11 were observed in the HTP imputed HLA dataset.
# Individuals with imputed heterodimer combination DQ2.5 + DQ2.2 + DQ7.5 were assumed to be carriers of HLA-DQ2.5 in trans, rather than true HLA-DQ2.5 (phasing error).
# Individuals with other combinations involving HLA-DQ2.5 were assumed to have true HLA-DQ2.5.
# Of course, these assumptions may be flawed given that this is imputed HLA data.

# 15 genotypes reported in the Sharp paper (* not observed in HTP):
# *DQ2.5/DQ2.5
# DQ2.5/DQ2.2
# DQ7.5/DQ2.2
# DQ2.5/DQ8
# DQ2.5/X
# DQ7.5/DQ2.5
# *DQ8/DQ8
# DQ2.2/DQ8
# *DQ2.2/DQ2.2
# DQ7.5/DQ8
# DQ8/X
# DQ2.2/X
# DQ7.5/X
# DQ7.5/DQ7.5
# X/X

df.HLA_DQ_Genotype <- df.Heterodimer_Combo %>%
  mutate(HLA_DQ_Genotype = ifelse(Heterodimer_Combo == "DQ2.5 + DQ2.5", "DQ2.5/DQ2.5",
                                  ifelse(Heterodimer_Combo == "DQ2.5 + DQ2.2", "DQ2.5/DQ2.2",
                                         ifelse(Heterodimer_Combo == "DQ2.5 + DQ2.2 + DQ7.5", "DQ7.5/DQ2.2", # This assumption has since been validated for the Celiac cases via molecular HLA typing.
                                                ifelse(Heterodimer_Combo == "DQ2.5 + DQ8", "DQ2.5/DQ8",
                                                       ifelse(Heterodimer_Combo == "DQ2.5 + X", "DQ2.5/X",
                                                              ifelse(Heterodimer_Combo == "DQ2.5 + DQ7.5", "DQ7.5/DQ2.5",
                                                                     ifelse(Heterodimer_Combo == "DQ8 + DQ8", "DQ8/DQ8",
                                                                            ifelse(Heterodimer_Combo == "DQ8 + DQ2.2", "DQ2.2/DQ8",
                                                                                   ifelse(Heterodimer_Combo == "DQ2.2 + .DQ2.2", "DQ2.2/DQ2.2",
                                                                                          ifelse(Heterodimer_Combo == "DQ8 + DQ7.5", "DQ7.5/DQ8",
                                                                                                 ifelse(Heterodimer_Combo == "DQ8 + X", "DQ8/X",
                                                                                                        ifelse(Heterodimer_Combo == "DQ2.2 + X", "DQ2.2/X",
                                                                                                               ifelse(Heterodimer_Combo == "DQ7.5 + X", "DQ7.5/X",
                                                                                                                      ifelse(Heterodimer_Combo == "DQ7.5 + DQ7.5", "DQ7.5/DQ7.5",
                                                                                                                             ifelse(Heterodimer_Combo == "X + X", "X/X",
                                                                                                                                    NA)))))))))))))))) %>%
  rename(Imputed_heterodimers = Heterodimer_Combo) %>%
  select(ID_colnames, EXCLUDE_from_analysis, EXCLUDE_reason,
         DQ2.5, DQ8, DQ2.2, DQ7.5, Imputed_heterodimers, HLA_DQ_Genotype)

#### Print the dataframe:
df.HLA_DQ_Genotype

Print the imputed heterodimers and their translation to an HLA-DQ genotype as reported in Sharp et al., 2019:

alleles_to_heterodimers


heterodimers_to_HLAgenotype

Output heterodimers_to_HLAgenotype to CSV and TSV

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

fwrite(heterodimers_to_HLAgenotype, "HLA_041822_Imputed_HLA_heterodimers_to_Inferred_HLADQ_Genotype_v0.1_JRS.csv")

Verify that each HTP_ID is only included in this dataframe once after EXCLUDE filter

df.HLA_DQ_Genotype %>%
  filter(EXCLUDE_from_analysis == 0) %>%
  group_by(RecordID) %>%
  summarise(N = n()) %>%
  arrange(desc(N)) %>%
  filter(N > 1)

print("Good.")

[1] "Good."

Verify that each participant has exactly one HLA-DQ genotype, and only one status per allele (0 or 1):

df.HLA_DQ_Genotype %>%
  select(MEGA.IID, HLA_DQ_Genotype) %>%
  group_by(MEGA.IID) %>%
  summarise(N = n()) %>%
  arrange(desc(N)) %>%
  filter(N != 1 | is.na(N))

df.HLA_DQ_Genotype %>%
  select(MEGA.IID, HLA_DQ_Genotype) %>%
  group_by(MEGA.IID, HLA_DQ_Genotype) %>%
  summarise(N = n()) %>%
  arrange(desc(N)) %>%
  filter(N != 1 | is.na(N))

`summarise()` has grouped output by 'MEGA.IID'. You can override using the `.groups` argument.

print("Good, each participant has exactly one HLA-DQ genotype, and only one status per allele (0 or 1).")

[1] "Good, each participant has exactly one HLA-DQ genotype, and only one status per allele (0 or 1)."

Prepare long-format dataset of HLA-DQ genotypes (1 if the person has the genotype, 0 if they have some other genotype)

df.HLA_DQ_Genotype.LONG <- df.HLA_DQ_Genotype %>%
  select(-c(DQ2.5, DQ8, DQ2.2, DQ7.5, Imputed_heterodimers)) %>%
  unique() %>%
  mutate(dummy = 1) %>%
  spread(key = HLA_DQ_Genotype, value = dummy) %>%
  gather(key = HLA_DQ_Genotype, value = Dosage, `DQ2.2/DQ8`:`X/X`) %>%
  mutate(Dosage = ifelse(is.na(Dosage), 0, Dosage))

df.HLA_DQ_Genotype.LONG

Prepare long-format dataset of HLA-DQ genotypes (1 if the person has the genotype, 0 if they have some other genotype)

temp1 <- df.HLA_DQ_Genotype %>%
  select(ID_colnames, EXCLUDE_from_analysis, EXCLUDE_reason, HLA_DQ_Genotype) %>%
  separate(HLA_DQ_Genotype, into = c("Heterodimer1", "Heterodimer2"), sep = "/", extra = "merge", remove = TRUE) %>%
  select(ID_colnames, EXCLUDE_from_analysis, EXCLUDE_reason, Heterodimer1) %>%
  unique() %>%
  mutate(Heterodimer_Index = 1) %>%
  rename(Heterodimer = Heterodimer1)

temp2 <- df.HLA_DQ_Genotype %>%
  select(ID_colnames, EXCLUDE_from_analysis, EXCLUDE_reason, HLA_DQ_Genotype) %>%
  separate(HLA_DQ_Genotype, into = c("Heterodimer1", "Heterodimer2"), sep = "/", extra = "merge", remove = TRUE) %>%
  select(ID_colnames, EXCLUDE_from_analysis, EXCLUDE_reason, Heterodimer2) %>%
  unique() %>%
  mutate(Heterodimer_Index = 2) %>%
  rename(Heterodimer = Heterodimer2)

df.HLA_heterodimers.LONG <- rbind(temp1, temp2) %>%
  group_by(FamilyID, RecordID, MEGA.FID,
           MEGA.IID, MEGA.LabID,
           EXCLUDE_from_analysis, EXCLUDE_reason, Heterodimer) %>%
  summarise(Dosage = n()) %>%
  spread(key = Heterodimer, value = Dosage) %>%
  gather(DQ2.2:X, key = Heterodimer, value = Dosage) %>%
  mutate(Dosage = ifelse(is.na(Dosage), 0, Dosage))

`summarise()` has grouped output by 'FamilyID', 'RecordID', 'MEGA.FID', 'MEGA.IID', 'MEGA.LabID', 'EXCLUDE_from_analysis', 'EXCLUDE_reason'. You can override using
the `.groups` argument.

df.HLA_heterodimers.LONG

Remove the exclusion columns; we will add them back once the analysis metadata is finalized

df.HLA_DQ_Genotype01 <- df.HLA_DQ_Genotype %>%
  select(-c(EXCLUDE_from_analysis, EXCLUDE_reason))

df.HLA_DQ_Genotype.LONG01 <- df.HLA_DQ_Genotype.LONG %>%
  select(-c(EXCLUDE_from_analysis, EXCLUDE_reason))

df.HLA_DQ_Genotype01

df.HLA_DQ_Genotype.LONG01

Output datasets to CSV and TSV

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

fwrite(df.HLA_DQ_Genotype01, "HLA_041822_HLA_Dosage_Heterodimers_DQgenotypes_SharpDQ7.5_v0.1_JRS.csv")
fwrite(df.HLA_DQ_Genotype01, "HLA_041822_HLA_Dosage_Heterodimers_DQgenotypes_SharpDQ7.5_v0.1_JRS.tsv", sep = "\t")

setwd(dir)
fwrite(df.HLA_DQ_Genotype.LONG01, "HLA_041822_HLA_Genotype_Dosage_LONG_v0.1_JRS.csv")
fwrite(df.HLA_DQ_Genotype.LONG01, "HLA_041822_HLA_Genotype_Dosage_LONG_v0.1_JRS.tsv", sep = "\t")

setwd(dir)
fwrite(df.HLA_DQ_Genotype.LONG01, "HLA_041822_HLA_Heterodimers_Dosage_LONG_v0.1_JRS.csv")
fwrite(df.HLA_DQ_Genotype.LONG01, "HLA_041822_HLA_Heterodimers_Dosage_LONG_v0.1_JRS.tsv", sep = "\t")